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Supersonic jet plumes were studied using a two-equation turbulence model employing corrections for com- 
pressible dissipation and pressure-dilatation. A space-marching procedure based on an upwind numerical scheme 
was used to solve the governing equations and turbulence transport equations. The computed results indicate 
that two-equation models employing corrections for compressible dissipation and pressure-dilatation yield im- 
proved agreement with the experimental data. In addition, the numerical study demonstrates that the computed 
results are sensitive to the effect of grid refinement and insensitive to the type of velocity profiles used at the 
inflow boundary for the cases considered in the present study. 


Nomenclature 

D = initial diameter of the jet 
k = turbulent kinetic energy 
M - Mach number 
P = static pressure 

P k ~ production of the turbulent kinetic energy 
R ~ gas constant 
r - radial coordinate 
T = static temperature 
/ = time 

u = velocity component in the axial direction 
x = axial coordinate 

a = model constants for compressibility and pressure- 
dilatation corrections in Eqs. (7) and (8) 
y = ratio of specific heats 
e — dissipation of turbulent kinetic energy 
fi = dynamic molecular viscosity 
p — density 

a - spreading rate of the jet, a = (dy, 2/ ,/dr) 
o) = vorticity 

Subscripts 
c = centerline 
j = jet 
/ = laminar 

t ~ turbulent 
* = freestream 

Introduction 

T HE recent resurgent interest in a High Speed Civil Trans- 
port (HSCT) and the National AeroSpace Plane (NASP) 
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clearly demonstrates the need for advanced propulsion sys- 
tems for supersonic velocities and beyond. Because of the 
complexity of the problem, numerous research programs have 
been initiated. One aspect of this research has been directed 
towards detailed understanding of the complex flowfield as- 
sociated with the engine over a wide range of operating con- 
ditions. Computational fluid dynamics (CFD) in conjunction 
with recent advances in turbulence modeling is being used 
extensively for detailed investigation of the propulsion related 
flowfield. Currently, there is a renewed interest in the flow 
properties of compressible mixing layer and jet plumes since 
it is central to many advanced propulsion systems. In addition, 
the classical fluid physics represented by this problem makes 
it ideally suited for basic numerical studies. 

An effort is currently underway to investigate the aero- 
dynamics of jet plumes, where the shear layer interacts strongly 
with the shock structure in much of the jet plume flowfield. 
It is widely accepted that the computational economy of the 
two-equation turbulence models offers the best compromise 
for computing complex turbulent shear flows. However, some 
deficiencies remain in these models that restrict their use in 
complex situations. In particular, the local density extensions 
of standard incompressible turbulence models were found to 
be inadequate in duplicating the experimentally observed re- 
duction in growth rate of the mixing layer with increasing 
convective Mach number. Various modifications’ 3 have been 
suggested, but none appear to model the physics of the flow 
adequately when parameters beyond mean growth rate are 
considered. Recently, Sarkar et al. 4 and Sarkar 5 recognized 
the importance of including compressible dissipation and pres- 
sure-dilatation effects in the two-equation turbulence models 
when computing high-speed flows. Simple corrections were 
proposed for compressible dissipation and pressure-dilatation 
that can be easily included in the existing two-equation tur- 
bulence models. In addition, Dash h included an axisymmetric 
correction in the turbulence transport equations when solving 
jet problems to account for vortex stretching. The standard 
two-equation model is recovered when the model constants 
for these corrections are assumed to be zero. 

During the past decade, considerable progress has been 
made in the area of CFD. Most of the research activities have 
been centered around developing efficient numerical algo- 
rithms for the solution of parabolized and full Navier-Stokes 
equations. Chief among these are upwind schemes because 
they are robust and efficient with good convergence prop- 
erties. Complex three-dimensional laminar and turbulent flows 
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over realistic aerodynamic configurations have been routinely 
simulated using these schemes. Despite the popularity of up- 
wind schemes, their extension to solving stiff equations are 
limited in literature, especially within the framework of space- 
marching procedures. 

Since the streamwise gradients are small in mixing layer/ 
jet plume problems, space-marching procedures are attractive 
alternatives to time-marching procedures. Lakshmanan and 
Abdol-Hamid 7 had studied the high-speed shear layer using 
an implicit upwind space-marching code containing the com- 
pressibility corrected two-equation turbulence model. The space 
marching was found to be highly efficient: two-orders of mag- 
nitude reduction of the residual was achieved in less than 15 
iterations in each marching plane using a CFL number as high 
as 15. 

The primary objective of the present work is to incorporate 
axisymmetric correction, compressible dissipation, and pres- 
sure-dilatation terms in the two-equation turbulence model. 
The resulting model is then utilized to compute high-speed 
compressible jet plumes. Another objective of the work is to 
investigate the effect of inflow profiles and grid refinement 
on the computed flowfield. 

Theoretical Formulation 

Governing Equations 

The theoretical formulation of the problem starts with the 
Favre-averaged form of the equations representing conser- 
vation of mass, momentum, energy, and turbulence quan- 
tities. The overbar is used to denote a conventional Reynolds 
average, whereas the tilde is used to denote the Favre-aver- 
age. For the sake of brevity, only the k and e equations are 
given here. Extensive details of the governing equations can 
be found in Ref. 8: 
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where P k = - pu \u\ (du,idx f ) is the production of the tur- 
bulent kinetic energy. Here, —pu' t u] is the Reynolds stress, 
and p k and /x, are defined as follows: 


simulation of isotropic turbulence. A simplified form of the 
axisymmetric correction used by Dash* is used in the present 
analysis. This additional term represents the contribution due 
to vortex stretching in axisymmetric flows. The details of the 
corrections are given as follows: 

Axisymmetric correction 



Model for compressible dissipation 

F. t = a l M;e s (7) 

Model for pressure-dilatation 

p"d" — — a 2 pP k M; -I- ot^pt\M; (8) 


where M, ~ V2 kiyRT is the turbulent Mach number, and 
e s is the solenoidal dissipation. Based on direct numerical 
simulation, Sarkar 5 recommends a, = 0.5, a, = 0.4, and 
a* = 0 . 2 . 

Boundary Conditions 

The boundary conditions to be imposed are problem-de- 
pendent. In the case of jet plume calculations, four types of 
boundaries are encountered. These include inflow, outflow, 
symmetry, and outer boundaries. Riemann invariants were 
used along the outer boundary corresponding to the external 
freestream. Symmetry boundary conditions were applied along 
the axis of symmetry. At the inflow boundary, three different 
initial velocity profiles were studied to investigate their effect 
on the computed flowfield. These include an experimentally 
measured velocity profile/' a hyperbolic tangent profile, and 
a top-hat profile. Since the solutions were obtained by a 
marching procedure, the downstream boundary conditions 
were left unspecified. 

For the turbulence transport equations either zeroth-order 
extrapolation or freestream values are used for k and e along 
the outer boundaries. If the flow is outgoing along the outer 
boundary, zeroth-order extrapolation is used. If there is flow 
entrainment, then freestream values are used along the outer 
boundaries. At the inflow laminar solutions were computed 
for the first three planes, and thereafter the k-e turbulence 
model is turned on. The inflow profile for k takes the same 
shape as the vorticity profile except it is multiplied by a spec- 
ified value of maximum turbulence intensity (Fig. 1). Once 
the k profile is known, the e profile is obtained using pro- 
duction equals dissipation hypothesis. Since the k-e equations 
were solved by a marching procedure, the downstream bound- 
ary conditions were left unspecified. 



p k = /x, + (pMk) ( 4a ) 

/A = Pt + (ft t lo- r ) (4b) 

P, = C^pikVe) (5) 

In Eq. (4), p, and p, represent the laminar and turbulent 
coefficients of viscosity, respectively. The model constants 
used in the analysis are C M = 0.09, C, = 1.44, C\ = 1 .92, 
C, = 0.79, o-, = 1.0, and = 1.3. In Eq. (1), and p”dT 
represent the contributions due to compressible dissipation 
and pressure-dilatation, and x is the axisymmetric correction. 
The compressible and pressure-dilatation effects are included 
using recently developed models 4 s based on direct numerical 




Fig. 1 Inflow boundary condition for k-e turbulence model: a) flat 
plate and b) jet flow . 
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Method of Solution 

The governing equations for the mean flow and turbulence 
quantities were integrated using an implicit upwind numerical 
procedure employing the approximate-factorization tech- 
nique. The numerical code used in this study has been de- 
signed to accommodate multiblock/multizone grid configu- 
rations. Two numerical algorithms are available in the code: 
1) van Leer flux-vector splitting and 2) Roe’s flux difference 
splitting. In addition, the code has options for three different 
algebraic turbulence models and a two-equation Jones-Laun- 
der turbulence model. 

Depending on the complexity of the problem, this code has 
the capability to compute the flowfield using either the space- 
marching or the time-marching procedures. It should be noted 
that when using a space-marching procedure for supersonic/ 
subsonic mixing problems, the upstream influence from all 
the points downstream is neglected. However, not all the 
streamwise derivatives are dropped. In particular, the stream- 
wise pressure gradient is fully retained in the supersonic re- 
gions of the flow. In the subsonic region of the flow, the 
streamwise pressure gradient term is dropped. This approach 
is valid as long as the pressure gradient is negligible. Extensive 
details of the numerical scheme, turbulence models, and the 
capabilities of the code can be found in Ref. 10. 

Results and Discussion 

Supersonic jet plumes (Figs. 2 and 3) exhausting into a 
nominally stationary external stream were computed using an 
improved form of the two-equation turbulence model. In all 
the cases, the initial jet radius was equal to 1 .28 cm. Because 
of the axisymmetric nature of the problem, only a sector of 
the jet was computed. The outer boundary of the computa- 
tional domain in the radial direction was located at eight 
radius from the jet axis. 

A compact algebraic grid generation program" was used 
for generating the mesh. The computations employed 401 
equally spaced points in the axial direction. A grid refinement 
study was carried out in the radial direction employing 81, 
101, and 121 points. Grid-independent results were obtained 
for computations employing 121 points. Approximately 61 
points were located in the initial jet radius with a typical grid 
spacing of 0.036 cm at the axis and 0.005 cm along the lip 
line. It should be noted that of the 61 points used in the jet, 
41 points were placed between 0- 1 . 18 cm, and the remaining 
21 points were placed in the narrow region between 1.18— 
1.28 cm. Using this grid resolution it was possible to inter- 
polate the velocity profiles available from the experimental 
measurements. g The remaining 61 points were distributed be- 



Fig. 2 Schematic of a fully expanded supersonic jet plume. 



Fig. 3 Schematic of an underexpanded supersonic jet plume. 


tween the lip line and the outer boundary using a cubic dis- 
tribution of grid cell sizes. 

The average computational time for each case is approxi- 
mately 140 CPU seconds on the Cray 2 computer. Approxi- 
mately 10-15 iterations were needed in each plane to achieve 
two orders-of-magnitude reduction of the residual using a 
CFL number as high as 15. The computed results did not 
differ significantly when the calculations were allowed to achieve 
three orders-of-magnitude reduction of the residual. 

For cases of an under or overexpanded jet flow, the cal- 
culation was carried up to 20 jet diameters in the streamwise 
direction. For fully expanded jet plumes the computation was 
carried up to 40 jet diameters in the axial direction. The 
freestream conditions for all the cases were 

M , = 0.05, P, = 101,325 N/m 2 , T , = 293 K 

Case I: Fully Expanded Supersonic Jet 

A fully expanded supersonic jet exhausting into the as- 
sumed freestream condition was computed. For this case, the 
jet operating conditions were 

M, = 2.2. PJP, = 1.0, TJT. = 1.0 

The computations were compared with the experimental 
measurements of Eggers. 1 ' Figure 4 shows the effect of axi- 
symmetric correction on the centerline velocity. These results 
are obtained using the turbulence model employing correc- 

With axisymmetric term [6] 1 

> With corrections [4,51 

Without axisymmetric term J 

O Experiment [9] 



Fig. 4 Effect of axisymmetric correction on the centerline velocity 
for a fully expanded supersonic jet plume. 


Inflow velocity profile 



Fig. 5 Effect of inflow profiles on the centerline velocity for a fully 
expanded supersonic jet plume. 
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P =101325 N/m 2 ,T = 293K,M =0.05 

oo J oo ^ oo 



Fig. 6 Effect of inflow profiles on the centerline turbulence intensity 
for a fully expanded supersonic jet plume. 


With corrections [4,5] 



X/D 

Fig. 7 Effect of compressibility and pressure-dilatation corrections 
on the centerline velocity for a fully expanded supersonic jet plume. 

P oo =101325 N/m 2 , T^r 2?3K, M = 0.05 



Fig. 8 Spreading rate for a fully expanded supersonic jet plume. 


tions for compressible dissipation and pressure-dilatation ef- 
fects. It is clear that by using the axisymmetric correction, 
slightly slower decay rate for the jet is observed. On the other 
hand, without the axisymmetric correction better agreement 
with the data for the velocity decay is obtained. In addition, 
numerical difficulties were experienced with the axisymmetric 
correction term when computing shock containing jet plumes. 
As a result, the axisymmetric correction term is neglected in 
subsequent computations. Figures 5 and 6 show the effect of 
inflow profiles on the centerline velocity and turbulence in- 
tensity. These results show that the centerline values are in- 
sensitive to the type of profiles used at the inflow boundary. 


Figure 7 presents the effect of compressibility and pressure- 
dilatation corrections on the computed centerline velocity. It 
is clear that the standard turbulence model provides a faster 
velocity decay, whereas the turbulence model with corrections 
yields a velocity decay that is in better agreement with the 
experimental data. The core length (distance required for 
mixing layer to reach the jet centerline) predicted by the 
standard turbulence model was about 11D, whereas the es- 
timated length using the model with corrections was about 
16D. Experimental results indicate a core length of about 
10Z). Figure 8 shows the jet spreading rate with and without 
compressibility correction. It is clear that the model without 
correction predicts a higher growth rate as compared to the 
compressibility corrected model. 

Figure 9 presents a comparison of the computed axial ve- 
locity profile at X/D = 13.5 and experimental data. The 
turbulence model with corrections predicts a velocity profile 
that is in good agreement with the experimental data, al- 
though the standard turbulence model significantly under- 
predicts the measured velocity profile in the center of the jet 
and overpredicts it near the edges. 

Figure 10 shows the normalized spreading rate for a fully 
developed flow (tr,, being the incompressible spreading rate 
equal to 0.1) as a function of the jet Mach number. It is clear 
that computations without any compressibility correction to 
the models (a = 0.0) show only a mild decrease in spreading 
rate with an increase in Mach number. However, inclusion 
of the compressibility correction to the models (a = 1 .0) lead 
to a sharp reduction in growth rate in the high Mach number 
regime. These results are in general agreement with Langley 
experimental data for compressible mixing layer. 


P_=101 325 N/m*. 1^= 293K, M oo = 0.05 



Fig. 9 Comparison of axial velocity profile at X/D = 13.5 with ex- 
perimental data for a fully expanded supersonic jet plume. 


With corrections [4, 5] 
- - - - Without corrections 



Fig. 10 Variation of the spreading rate with the jet exit Mach num 
ber. 
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Case 2: Underexpanded Supersonic Jet 

An underexpanded supersonic jet was computed for an 
identical freestream condition as in the previous case: In this 
case, the jet operating conditions were 

M, = 2.0, PJP, ~ 1.4, T JT, - 1.0 

Figure 11 shows the computed results using the turbulence 
model and comparison with the experimental data. 1 ■■ The com- 
puted results are shown for the first 20 jet diameters in which 
8 shock cells occur. The computations slightly underpredict 
the peak pressure in the first four cells of the jet. However, 
the shock cell spacing is predicted well. Figure 12 shows the 
effect of turbulence model with and without corrections on 
the computed pressure distribution along the jet axis. The 
results obtained using the standard turbulence model and the 
turbulence model with corrections are in good agreement for 
the first three cells. Beyond these three cells, the standard 
turbulence model predicts a faster decay of the remaining 
shock cell structure. Figure 1 1 also shows the effect of grid 
refinement on the computed centerline pressure distribution. 
It is clear that results obtained using refined mesh are in better 
agreement with the experimental data than the coarse mesh 
results in which the jet shock structure is displaced slightly 
towards the right. 

The computed rms values of the streamwise turbulence 
intensity shown in Fig. 13 along the jet lipline indicate a large 
improvement in the results when compressibility and pres- 
sure-dilatation corrections are included. Here, the experi- 
mental data are available only up to 10 jet diameters. In 
general, the agreement of computed results using the standard 
turbulence model with the experimental data was not good. 
The measurements and the computed results indicate that the 
turbulence intensity amplifies to reach peak values at the end 

Grid 

401 x 121 I 

\ With corrections [4,5] 

401 x 81 J 



Fig. 11 Comparison of centerline pressure distribution for an un- 
derexpanded supersonic jet plume. 



Fig. 12 Effect of compressibility and pressure-dilatation corrections 
on the computed centerline pressure distribution for an underex- 
panded supersonic jet plume. 


With corrections [4,5] 
Without corrections 



Fig, 13 Effect of compressibility and pressure-dilatation corrections 
on the lip line turbulence intensity for an underexpanded supersonic 
jet plume. 


With corrections [4,5] 



Fig. 14 Comparison of centerline pressure distribution for an over- 
expanded supersonic jet plume. 


With corrections [4,5] 



Eig. 15 Comparison of pressure distribution along the radial line 
rll) - 0.375 for an overexpanded supersonic jet plume. 


of compression zones and minimum values at the end of ex- 
pansion zones. 

Case 3: Overexpanded Supersonic Jet 

An overexpanded supersonic jet was computed for identical 
conditions given in the previous case except for the difference 
in pressure ratio at the nozzle exit. In this case, the jet op- 
erating conditions were 


M , = 2 . 0 , 


rjr - l.o 


PJP, - 0.7, 
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The computed results for the pressure distribution along 
the jet axis and a radial line rlD = 0.375 are shown in Figs. 
14 and 15, respectively. The centerline pressure distribution 
shown in Fig. 14 represents the axial extent of available data 
for this case. The numerical solution predicts the irregular 
looking pressure data 1 - 14 reasonably well. The comparison of 
the computed pressure distribution and the experimental data 
along the outer radial line is very good for the first three 
shock cells. The location and magnitude of the peak pressure 
were predicted well in the computations. In the downstream 
region, significant differences were observed in the shock wave 
length between the computed results and the experimental 
values, suggesting that measured jet plume has a faster jet 
spreading rate. 

Concluding Remarks 

Compressible high-speed jet plumes were studied using a 
two-equation turbulence model. It was observed that for 
supersonic/subsonic mixing problems, converged results could 
be obtained efficiently by using the space-marching proce- 
dure. In addition, corrections for compressible dissipation and 
pressure-dilatation for the two-equation turbulence model were 
tested for high-speed jet plumes. For supersonic jets, the 
turbulence model with corrections had significantly improved 
the agreement between computed results and experimental 
data. Moreover, the numerical study demonstrated that the 
computed jet decay rate was sensitive to grid resolution and 
insensitive to the type of profiles used at the inflow boundary. 
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